Decay properties of the density matrix and Wannier functions for interacting systems 
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For non-interacting electrons the one-particle density matrix and the related Wannier functions 
characterize a material as insulating or metallic. Introducing many-body Wannier functions, we 
show that this characterization can be carried over to interacting systems. In particular, we find 
that an exponential decay of the density matrix characterizes not only band insulators but also Mott 
insulators. The properties of the many-body Wannier functions differ, however, from those of the 
Wannier functions of a non-interacting systems. 
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The locality properties of solids within independent- 
particle theories have been at the focus of much attention 
It was shown that for insulators the density matrix 
decays exponentially. The rate of the decay is related 
to the size of the energy gap in the band-structure: A 
larger gaps mean a faster decay. In contrast, for a metal 
at zero temperature, the decay is much slower, namely 
algebraic. The decay properties of the density matrix are 
reflected in the locality properties of the Wannier func- 
tions [|||^ . The Wannier function Wj (r — R) associated 
with the j-th band is localized around the lattice vector 
R, i.e. it tends to zero for r far from R. In an insula- 
tor this decay is exponential for isolated bands Q]. The 
practical importance of the Wannier functions is two- 
fold. First, they allow for an intuitive interpretation of 
polarization effects and bonding properties in a solid %. 
Localized Wannier functions correspond to either bonds 
or lone electron pairs. Second, they are the basic quan- 
tity in some linear-scaling algorithms |6[0]. When they 
are localized, they need only be calculated within the lo- 
calization region, allowing thus to use divide and conquer 
strategies that finally lead to linear scaling. 

In the following we will show how the characterization 
of an insulator by the exponential decay of the density 
matrix can be carried over to interacting systems. Using 
natural orbitals and what we call many-body Wannier 
functions allows us to discuss these many-body proper- 
ties in terms of single particle wave-functions. It turns 
out that the localization properties of the occupied many- 
body Wannier functions are similar to those of the Wan- 
nier functions for non-interacting systems. For the vir- 
tual Wannier functions there are, however, qualitative 
differences. 

Criteria for characterizing the insulating state that are 
valid beyond the independent-particle picture have al- 
ready been put forward. Kohn has shown that the 
many-body wavefunction of an insulating ring breaks up 
into a sum of functions which are localized in discon- 
nected regions of the many-particle configuration space. 
Resta and Sorella |^ have shown that the expectation 
value of the operator ^^i"^"^!^)^ can be used to distinguish 



between insulators and metals. This criterion seems, 
however, to be restricted to one-dimcnsional systems. As 
we will show, a more general criterion is given by the 
decay of the density matrix. All these criteria have in 
common that they refer only to the ground state wave- 
function and do not require information about excited 
states. Criteria of this type are the most useful ones in 
numerical calculations, where in general only the ground 
state is calculated. 

We start our construction by generalizing the concept 
of Wannier functions to interacting systems. Let us con- 
sider a iV-electron wavefunction ^'(xi,X2, ...,XAr), where 
X — {r, s} is a combined spatial and spin coordinate. The 
reduced density matrix -D(r, r') is then given by 



D(r,r') 



N \ dsdx'}...dx 



N 



**(r, s, X2, ...xjv) *(r', s, X2, ...xat) 



(1) 



The eigenfunctions $i(r) of the reduced density matrix 
are called the natural orbitals, its eigenvalues rii the nat- 
ural occupation numbers. In this basis the many-body 
density matrix takes the same form as an independent- 
particle density matrix: 

i 

Translating all the spatial arguments of '5 by a lattice 
vector R multiplies the wavefunction by a phase factor. 
Since these phase factors cancel by the definition of the 
density matrix D we find 

L'(r + R,r' + R) = D{r,r') . 

D{r',r) has thus the same symmetry properties as 
the Hamiltonian of a periodic solid in an independent- 
particle picture Therefore its eigenfunctions can be 
written as Bloch functions, labeled by a band index j and 
a vector k contained in the Brillouin zone 



(2) 



The Uj^k{r) are periodic functions with respect to the real 
space primitive cell. From these functions we construct 
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the k-dependent density matrix, which plays a central 
role in our argument: 

Dk(r,r')=^n,(k)C/,-k(r)!7;_k(r'), 

J 

where the rij (k) are the occupation numbers correspond- 
ing to the natural orbitals (Q) . The full density matrix is 
related to the k-dependent density matrix by 

I?(r, r') = / rfk D^{v, r') e^'^^"-""') , (3) 

(27r)-^ Jbz 

where V is the volume of the real-space primitive cell 
and the integration is over the Brillouin Zone (BZ). By 
construction the J7j_k(r) satisfy the eigenvalue equation 



j dr' Dkir,r') Uj.y,{v') - nj(k) t/,-u(r) 



(4) 



Because of the analogy to the band structure problem, we 
call the eigenvalues rijik) the occupation band-structure. 
It has the periodicity of the Brillouin zone, and N- 
representability requires that nj{k) G [0,2] ([0,1] 
without spin degeneracy). 

We next use the natural orbitals (||) to define many- 
body Wannier functions, Wj (r — R) 

W,{v _ R) = J:_ /■ dk e'''(--^) C/,,k(r) . (5) 

l^^^J JBZ 

This definition is formally similar to the usual definition 
of the independent-particle Wannier functions: 



K'^T^) JBZ 



where (^j.K(r) are the independent-particle Bloch func- 
tions 

0j,K(r) = e**'''wj-k(r), 
whose cell-periodic part satisfies the eigenvalue problem 



j dr' i?0(r, r') u,- K(r') = e, (k) w,- k(r) 



(6) 



We notice that in the non-interacting case the eigenvalue 
problem (Q) involving the density matrix is replaced by 
(^, involving the effective Hamiltonian _ff^(r, r'). The 
reason for this is that in the limit of no interaction the 
density matrix becomes highly rank deficient, with zero 
occupation number for all virtual orbitals and all the oc- 
cupied orbitals being degenerate. This excludes, for the 
independent-particle Wannier functions, the use of the 
mathematical theorems onto which the following discus- 
sion is based, and also implies that in general the func- 
tions obtained from the many-body Wannier functions 
in the limit of no interaction differ from the usual non- 
interacting Wannier functions. 



Also the expression for the many-body density matrix 
in terms of the many-body Wannier functions is slightly 
more complicated than the corresponding expression in 
the independent-particle case. 

D{r,r')^ nj(R--R.')W*{r -Ii)Wj{r' -K') , 

where nj(R-R') = V/{2Tr)^ J dke'^^^-'^">nj(k), while 
for independent particles only the term nj(0) survives. 

From a mathematical point of view we can now dis- 
tinguish two cases. In the first case _Dk(r,r') is an an- 
alytic function of k, in the second case Dk(r, r') has a 
discontinuity at some k = kp. Then, in the absence 
of degeneracies, in the first case, the occupation band- 
structure rij (k) will be an analytic in k , while in the 
second case it will have a discontinuity at kp. We asso- 
ciate the first case with an insulator, while we call sys- 
tems with a discontinuity at kp a metal. To justify this 
classification we remark that experimentally the sharp- 
ness of the Fermi surface in metals is strongly suggested 
by the de Haas-van Alphen effect, cyclotron resonances, 
and Friedel oscillations. Additional support comes from 
the momentum distribution N{p). It can be probed by 
Compton scattering experiments [0. If N{p) is smooth 
the system is an insulator, while for a metal it shows 
a discontinuity. For jellium this discontinuity has been 
thoroughly investigated by many-body techniques [p^ . 
Being a one-particle property, the momentum distribu- 
tion can be determined from the reduced density matrix 

iV(p) = J drdr' e-'P^"-"") D{r,r') . 

Writing the density matrix in terms of the natural or- 
bitals one obtains 



iV(p)=^n,(p-G) 



p-G 



(G) 



where a^(G) are the plane wave expansion coefficients 
of [/j k(r). Obviously A^(p) is analytic if rijCk) and 
are analytic, whereas a discontinuity in nj(k.) manifests 
itself in a discontinuity in N{p). 

We are now in the position of making statements about 
the decay properties of the density matrix. For an insu- 
lator, where the density matrix is analytic in k, we find 
from (^) that, since the Fourier transform of an ana- 
lytic function decays exponentially [p^ , the density ma- 
trix D{r,r') decays exponentially fast to zero with in- 
creasing distance between r and r'. For a metal, on the 
other hand, Dk(r,r') has a discontinuity on the Fermi 
surface and hence the density matrix only decays alge- 
braically. The decay properties of the density matrix for 
interacting systems are thus qualitatively the same as for 
non-interacting systems. We emphasize that this results 
is valid for zero temperature, while at finite T the decay 
should become exponential also for a metal [fi5[. 
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We next turn to the properties of the Wannier func- 
tions using the same mathematical theorems. For an in- 
sulator, the eigenfunctions JTj.k are analytic functions of 
k, inheriting their analyticity from I?k(r, r') by cqn. (|^). 
Thus by eqn. (||) all the many-body Wannier functions 
will be exponentially localized for nondegenerate occu- 
pation bands. To make more specific predictions about 
the qualitative shape of the Wannier functions and to 
compare to the independent-particle case we will for the 
moment specialize our discussion to weakly correlated 
insulators, i.e. correlated insulators derived from band 
insulators. Then the occupation numbers of the (for- 
merly) occupied natural orbitals are close to two (one 
without spin degeneracy), while the occupation numbers 
of the (formerly) virtual natural orbitals are close to zero. 
Hence the occupation band structure is fairly flat. For 
a one-dimensional band insulator the occupation band 
structure has the qualitative form shown in Fig. |l|. 
2 
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FIG. 1. Qualitative sketch of the occupation band struc- 
ture Tijik) for a weakly correlated one-dimensional band in- 
sulator with lattice constant a. 



We will now put forward some qualitative predic- 
tions about the shape of the many-body Wannier func- 
tions of weakly correlated band insulators. As in the 
independent-particle case such a discussion is made dif- 
ficult by the fact the the Wannier functions are not 
uniquely defined. Multiplying the natural orbitals by a k 
dependent phase <I>j,k(r) $j,k(r) exp(i6'(k)) generates 
a new set of valid natural orbitals leading via eqn. (|^) to 
a new set of Wannier functions. In the following discus- 
sion we will assume that the phase 6{k) is chosen such as 
to minimize the overall variation with respect to k, i.e. 
we consider maximally localized Wannier functions jsj. 

Eqn. (^) shows that the variation of nj(k) is coupled 
to the one of <I'j^k(r)- We use here the word variation 
in a loose sense, meaning that either the function or 
its derivatives are large. Alternatively, a function has 
a strong variation if there are regions in which the con- 
vergence radius of its Taylor expansion is small. The 
fact that the occupation band structure is smooth con- 



sequently implies that the k variation of $j^k(r) is slow 
as well. The typical length scale over which the variation 
takes place in k space is ir/a. The decay constant of the 
many-body Wannier functions is therefore of the order of 
the interatomic distance independent of the index j. 

The occupied independent-particle Wannier functions 
have the same qualitative decay behavior. This follows 
by considering the eigenvalue equation (^ . Since the oc- 
cupied bands have in general a slow variation over the 
whole Brillouin zone, the decay constant is again of the 
order of the interatomic distance. This is also what is 
actually seen in numerical calculations showing that the 
occupied natural orbitals <I>j_k(r) are very similar to the 
occupied Kohn-Sham functions (/>j,k(r) |0. 

The virtual independent-particle Wannier functions, 
however, are more extended than the occupied 
independent-particle Wannier functions. This follows 
from the fact that the higher virtual bands tj (k) become 
free particle like and have a strong variation close to the 
boundaries of the Brillouin zone. This trend was also 
demonstrated by explicit calculations . 

We conclude that the occupied independent-particle 
and many-body Wannier functions of weakly correlated 
band insulators are very similar. The virtual Wannier 
functions, however, differ in the two contexts. The many- 
body Wannier functions are all localized within the same 
volume, whereas their independent-particle counterparts 
become more and more extended with increasing energy. 
A similar effect is well known in atoms, where the natural 
orbitals are all localized within the same volume whereas 
the virtual independent-particle orbitals form a Rydberg 
series with larger and larger extent. 

As an example of a strongly correlated insulator we 
consider a one-dimensional Hubbard model [|l8| at half 
filling. This system is a Mott insulator [|9| for any finite 
value of the on-site repulsion U . We have calculated 
the density matrix for a finite system of 12 sites with 
periodic boundary conditions. The results are shown in 
Fig. 1^. The exponential decay of the density matrix is 
clearly visible, being faster the stronger the correlation, 
i.e. the larger the Hubbard term U . The size of the gap 
for the systems shown in the figure are given in table ^. 

The occupation band-structure associated with the two 
extreme U -values of Fig. |^ is shown in Fig. ^. We notice 
that, again, a strong variation of the occupation band- 
structure implies a slower decay of the density matrix. 
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FIG. 2. Exponential decay of the density matrix Dij as a 
function of the inter-site distance j — i (in units of the lattice 
constant) for a Mott insulator. 



We have shown that the characterization of metals 
and insulators by the decay properties of the density 
matrix can be carried over to interacting systems. In 
addition we have defined many-body Wannier functions 
and discussed their localization properties. We expect 
that such Wannier functions will turn out to be useful 
in the construction of many-body 0(N) schemes. Lo- 
calization concepts arc a basic ingredient for the under- 
standing of bonding in solids and molecules. Our work 
extends these concepts that are well established in non- 
interacting framework to true many-body systems. 

S.G. thanks Nick Trefethen for pointing out references, 
and T. Arias, P. Horsch, O. Gunnarson, G. StoUhoff, J. 
Hutter, and D. Vanderbilt for interesting discussions. 
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U = exp(n/2) 
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4.4817 


2.090 


5 


12.1825 


8.840 


7 


33.1155 


29.385 


9 


90.0171 


86.142 



TABLE I. Gap Eg = E{N - 1) + 2E{N) - E{N + 1) for 
a periodic Hubbard chain of 12 sites for various values of U. 
All energies are given in units of the hopping matrix element, 
i.e. the band-width of the non-interacting system is 4. 
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FIG. 3. Occupation band structure of Mott insulators cor- 
responding to two different values of U. Clearly the variation 
of n{k) is stronger for smaller values of U. This is true even 
if the two curves are scaled such that they coincide at the 
end points. For smaller correlation n{k) is flat at the edges 
and decays mainly around kp , over a distance that is roughly 
half of the width of the cell, while for larger U the decay is 
nearly uniformly over the whole cell. The shorter length scale 
of the first curve manifests itself in a larger decay length of 
the density matrix. 
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